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ABSTRACT 

An analytic relation between the statistics of photons in pixels and the number 
counts of multi-photon point sources is used to constrain the distribution of gamma- 
ray point sources below the Fermi detection limit at energies above 1 GeV and at 
latitudes below and above 30°. The derived source-count distribution is consistent 
with the distribution found by the Fermi collaboration based on the first Fermi point 
source catalogue. In particular, we find that the contribution of resolved and unre- 
solved active galactic nuclei (AGN) to the total gamma-ray flux is below 20 to 25%. 
In the best fit model, the AGN-like point source fraction is 17 ± 2%. Using the 
fact that the Galactic emission varies across the sky while the extra-galactic diffuse 
emission is isotropic, we put a lower limit of 51% on Galactic diffuse emission and 
an upper limit of 32% on the contribution from extra-galactic weak sources, such as 
star-forming galaxies. Possible systematic uncertainties are discussed. 
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Introduction 



Faint gamma-ray point sources cannot be detected individually, however their presence af- 
fects the statistics of photons across the sky. One can use the observed statistics of photons to 
infer some general properties about the population of point sources below the Fermi detection 



sources is familiar in radio and X-ray obseryations (e.g., 



1974 



Hasinger et al. 



1993 



Mivaii & Griffiths 



2002 



Soltan 



Scheuer 



1957 



the stuc 


y of faint point 


Condon 


1974 




Scheuer 



201 ih . A closely related subject is the 



use of extreme stati s tics to constrain the non-Gaussianity of cosmic microwave background data, 



Colombi et al 



(hoil ) and references therein. In this paper, we extend the prior analysis 
( ISlatyer fc Finkbeinerl " 2010 ) and use the points in cells statistics to understand the population 



of gamma-ray point sources. 

In a larger context, the problem is to separ ate different source s of gamma-ray emission. 
The standard strategy is either to use templates ( IDobler et al.ll2010f ) that trace the sources or 
else to simulate the cosm ic ray propagation and gamma-ray production using, e.g., Galprop 



( jStrong et al.l 120071 . [2009J). In this paper we would like to adopt a different methodology: we 
assume some general properties of the sources in order to obtain model-independent constraints 
on the contribution from these sources. 

We separate the sources based on their statistical properties. In the paper we consider three 
sources of gamma-rays at high latitudes. The first source is a diffuse source that varies over 
large angles. The main contribution to this source comes from the Galactic diffuse emission 
(ir° production, ICS photons, and bremsst raMung) and, possibly, a population of faint Galactic 



point sources, such as millisecond pulsars (IFaucher-Giguere fc Loeb 



201 



0; 



Malyshev et al. 



2010 



?). We will call this source non- isotropic Galactic diffuse emission. It puts a lower limit on the 
actual Galactic emission. 

The second source corresponds to an isotropic distribution of gamma-rays, i.e., the statistics 
of photons across the sky is consistent with the Poisson distribution for this source. The isotropic 
flux has contributions from the homogeneous p art of the Galactic diffus e emission, from diffuse 
intergalactic emission (references can be found in lStecker fc Ventersll2010l ). and from a population 
of very weak extragalactic sources that on average emit much less than one photon during the 



Pavlidou fc Fields 



2002 



Fields et al. 



time of observation, such as star-forming galax ies (e.g. 

20101 ). starburst galaxies ([Thompson et al.l 120071 ) . galaxy clusters (jBerrington fc Dermerll2003l ). 
etc. The value of the isotropic emission puts an upper limit on the gamma-ray flux from very 
weak extra-galactic sources, as well as an upper limit on the contamination from cosmic rays 
(Ubdo et al.ll2010li l. 



-3- 



The third source is a population of point sources modeled by a broken power-law source- 
count distribution. We assume that the point sources are distributed homogeneously over the 
sky. The statistics of photons coming from these point sources has a non-trivial form (derived in 



Appendix El) different from 


re Poisson statistics. 


These sources mode 


a po 


Dulation of AGN-like 


point 


sources ( 


Stecker et al. 


1993; 


Padovani et al 


1993; 


Chiang et al. 


1995 


Stecker & Salamon 


1996 


Miicke & Pohl 


2000 


; lAbazajian et al. 


2010; 


Neronov & Semikoz 


2011) 





We note that in our approach it is impossible to separate the isotropic part of Galactic emis- 
sion from the diffuse extragalactic background (EGB). Consequently, all results will be quoted 
either with respect to the total gamma-ray flux or in absolute values. We find that at high lati- 
tudes (below and above 30°) for energies above 1 GeV the contribution of AGN-like point sources 
(both resolved and unresolved) to the total gamma-ray flux is 17% ± 2%, the contribution from 
Galactic diffuse emission is above 51% and the contribution from extra-galactic weak sources 
is below 32%. Using these values , we estimate that the contribution of unresolved AGN-like 
point sources to the EGB model in Abdo et "all ( 2010b ) is below ~ 25%, which is consistent with 



Abdo et al. 



fl2010d ). 



The paper is organized as follows. In Section [2} we describe a general model and a fitting 
algorithm. In Section [31 we present an analysis of the Fermi data. Section H] has discussion. In 
Appendix [A] we derive the statistics of photons coming from a population of point sources. In 
Appendix [HI we determine a model for non-isotropic Galactic diffuse emission. In Appendix [Cj 
we repeat the data analysis of Section [3] for different pixel sizes. 



2. Model 

In this section we describe the model that we later use to fit the Fermi gamma-ray data. At 
first, we present an analytic relation between the counts of photons in pixels and the statistics of 
point sources. Then we take into account the detector point spread function (PSF) and describe 
a model of Galactic diffuse gamma-ray emission. The model that we present in this section is 
rather general. An application of this analysis for the Fermi data is presented in Section [3j 



2.1. Statistics of Photon Counts and Point Sources 

We will assume a pixelation of the sphere with pixels of equal area. Usually a model of 
gamma-ray emission is represented in terms of the expected number x p of gamma-rays in pixel 
p. The problem is that the number and the positions of weak point sources are not known. The 
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position of a point source can be found only by a "detection" of the source. However, if one is 
interested only in the total number of sources of a given flux, then a "detection" of sources is 
not necessary. In this paper, we use the statistics of photon counts in pixels, in order to find the 
source-count distribution of point sources without actual identification of the sources. 

In order to find the statistics of photon counts in pixels we calculate the total number of 
pixels rik that contain k photons. In this calculation, the information about the position of the 
pixels on the sky is not preserved, but the remaining information may be sufficient to constrain 
general properties of the population of sources. There are two competing conditions for this 
method to work. On the one hand, there should be sufficiently many pixels to determine the 
statistics of photon counts, while, on the other hand, the size of the pixels cannot be too small 
compared to the PSF, otherwise the statistics corresponding to the presence of point sources 
cannot be distinguished from non-isotropic diffuse emission. 

If the total number of pixels is iVpi x and is the observed number of pixels with k photons, 
then we can estimate the probability to find k photons in a pixel as 

Pk = (1) 

1 v pix 

For large -/V P i x , the statistical uncertainty of n k is approximately y/nk (there is a small correction 
due to the fact that the total number of pixels is fixed, i.e., the process is multinomial). If 
there are no point sources and all the emission comes from an isotropic diffuse source, then 
the probability distribution pk is the Poisson probability of getting k photons given the mean 
rate. In the presence of point sources the probability distribution is more complicated than the 
Poisson distribution. In the context of X-ray point sources, the probability distribution is called 
a P{D) distributi on (or a P(D) diagram) and is usually computed with the help of Monte Carlo 



simulations (e.g.. iMiyaji fc Griffiths! 120021 ) . 



In this paper we find an analytic relation between the probability distribution of photon 
counts in pixels and the source-count distribution. We will use this relation to find a source- 
count distribution that has the best fit to the observed probability distribution determined in 
Equation ([I]). 

In calculation s we use the met hod of probability generating functions (an introduction can 



be found in, e.g., Ilioel et al.l [l97l| , Section 3.6). For a given discrete probability distribution 



Pk, k = 0,1,2..., the corresponding generating function is defined as a power series in an 
auxiliary variable t 

oo 

P(t) = J2p k t k - (2) 



k=0 
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If the generating function P(t) is known, then the probabilities pk can be found by picking the 
coefficient in front of t k , or, equivalently, by differentiating with respect to t 



Pk 



1 d k P{t) 



k\ dt k 



(3) 



t=o 



An important property of probability generating functions is that a sum of two independent 
random var i ables corresponds to a product of the corresponding probability generating functions 
(IHoel et al.l Il97ll ). For example, if there are two independent sources of gamma-rays on the 
sky, e.g., Galactic and extragalactic, then the probability generating function for the photon 
counts with the two sources overlaid is given by the product of probability generating functions 
corresponding to the two sources (see Appendix El for more details). 

Let x m denote the average number of sources inside a pixel that emit exactly m photons 
durin g the time of obs ervation. Using the sum-product property of probability generating func- 
tions (IHoel et al.lll97ll Section 3.6), we derive in Appendix |A] the following relation between the 



probability generating function for the photon counts and the expected number x m of m-photon 
sources 

oo / oo \ 

J2Pkt k = exp I ^{x m t m -x m )\. (4) 

k=0 \m=l J 

This formula provides an analytic relation between the expected numbers of m-photon sources 
and the statistics of photons in pixels. The probabilities pt to observe k photons in a pixel are 
determined by expanding the right hand side of this equation and picking the coefficient in front 
of t k . 

If we substitute t = e lU) , then Equation (121) becomes a disc rete Fourier transform analog of 



the probability characteristic functions (e.g.. IHoel et al.lll97ll ). The probability characteristic 
functions wer e extensively used in the study of radio point sources below detection limit (e.g., 



Scheuer 



19571 ). Similarly to probability generating functions, the characteristic function for a 



sum of two independent random variables corresponds to a product of the two characteristic 
functions. 

The statistics of point sources can be described in terms of the differential source counts 
as a function of flux S, denoted by dN/dS. One of the simplest forms for the source-count 
distribution is a broken power-law 




S > Sbreak 
S < Sbreak; 



(5) 
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where we have to require n\ < 2 and n 2 > 2 in order to have a finite total fluxo 

Let S be an average flux from a source (normalized to give the number of photons during 
the time of observation). Then the probability to observe m photons from this source is 

Cm 

Pm(S) = ^je~ s . (6) 
ml 

If the sources are distributed isotropically, then the average number of m-photon sources inside 
a pixel is given by the Poisson probabilities for sources with flux S to emit m photons 

x m = — dS—{S)—e , 7 
Air J db ml 

where f2 pix is the angular area of the pixel. 



2.2. PSF 

In the case of non-zero PSF a source in some pixel may contribute gamma-rays to nearby 
pixels as well. In order to correctly estimate the number of point sources of a certain strength, 
we need to know how often a gamma-ray from a source in one pixel is detected in a different 
pixel. 

We can represent the effect of the PSF as a smearing of a point source over some area, so 
that the flux from the point source is split between several pixels. We determine the average 
properties of the flux splitting from Monte Carlo simulations. In the simulation, we assume that 
the profile of the PSF is Gaussian. In order to calculate the average splitting, we put n Gaussians 
at random positions on the sphere and integrate every Gaussian over the pixels to get a collection 
of fractions /j, i = 1, . . . , N pix , such that f\ + f 2 H — • = 1 (in practice, one can truncate at some 
minimal value of /). Denote by An(f) the number of fractions for n Gaussians that fall within 
Af, then the average distribution of fractions is defined as 



p(/) - A " (/) 



nAf 



(8) 

A/— »0, n— >oo 



1 A benchmark model is provided by a static universe filled homogeneously with sources that have intrinsic 
luminosity Lq. The number of sources with flux larger than S — is N(S' > S) ~ R 3 ~ S~ 15 . The differential 
number count is dN/dS ~ —S~ 2 ' 5 . Index n = 2.5 provides a benchmark value for differential source count at 
large S. A break or a cutoff at small flux (large distance to the source) is necessary given the finite size of the 
visible universe. 
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This distribution is normalized as 

f fp(f)df = 1. (9) 



The case of zero PSF corresponds to p(f) = 8(f — 1). Sums over fractions are substituted by the 
integral f p(f)df. In particular, the contribution of point sources with intrinsic flux between S 
and S + dS to the total number of m-photon sources is 

dx m = dN(S) [ dfp{f)^^e-f s . (10) 
J ml 

The expected number of m-photon sources inside a pixel including the effect of PSF is 

x m = ^f dS^(S) [ df P (f)^e^ s . (11) 



An J Q dS J m\ 



An example of the function p(f) relevant for the data analysis in Section [3] is presented in Figure 
[TJ If we rescale the integration variable S in Equation (TTTT) by /, then the effect of the PSF can 
be represented as a change of the source-count function 

§(^>^,ff (,//). 

The apparent number of sources with small fluxes in ^ is larger than the corresponding number 
of physical sources in ^ due to contribution from point sources in nearby pixels, whereas the 
apparent number of sources with large fluxes is smaller due to loss of the flux to nearby pixels. 
One can see the effect of PSF in Figure [2] on the right. On the same figure, we also plot the 
expected number of m-photon sources for discrete m, derived from 



2.3. The Model 

In Equation (j4j), the probabilities pk and the average expected number of m-photon sources 
may depend on the position of the pixel. The quantity that we use is the averaged probability 
p k = -rr*-, where n k is the number of pixels with k photons. The generating function for averaged 
probabilities is given by 

$>^ fc = at- E exp 5>™ tm - o ■ ( 13 ) 

k=0 iVpix p=l \m=l / 

This form of the generating function can be used when the distribution of sources is not isotropic 
or the exposure is non-uniform. If the exposure is sufficiently uniform, then we may assume that 
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the distribution of extragalactic sources is independent of the pixel position and can be taken 
out of the sum over pixels. The diffuse emission corresponds to "one-photon" source density that 
depends on the pixel, x^mip)- Thus, the expression on the right-hand side of Equation f JT3|) can 
be split into a product 

^ Q X ditt t ~ x ditt+^m=l( Xmtm ~ Xm ) — J ^ g z diff f ~ x diff I . (iJ2m=l( x mt rn - x m) (14) 

N p [x P =i \ n p' ix P =i y 

The second term in the product on the right hand side corresponds to the isotropic distribution 
of AGN-like point sources 

P{t) = exp ( J2 (x m t m - x m ) ) , (15) 

\m=l J 

while the first term corresponds to diffuse emission or, indistinguishably, one-photon sources 

N ■ 

-i lv pix 

D{t) = e^*"^. (16) 

P =i 

It is convenient to separate the diffuse emission into a non-isotropic part that puts a lower limit 
on Galactic diffuse emission and an isotropic part that consists of isotropic Galactic component 
and a possible population of weak extragalactic sources (additional to AGN-like sources) 

^diff = X Gal + ^isotr- (17) 

The main problem with the non-isotropic diffuse emission is a variability on large scales (on small 
scales the variability in Galactic emission is subdominant to Poisson noise and variability due 
to the presence of undetected point sources). In Appendix [B] we construct a model for the non- 
isotropic component Xq^ that varies on large scales. We mask the known point sources and use 
low multipole spherical harmonics in order to filter out the small-scale variations while preserving 
the large-scale ones. We add a constant to x^ al so that min(xQ al ) = 0. In the following, the 
non-isotropic component of diffuse emission is fixed. The corresponding generating function is 

G(t) = —J2 e^al^Ga!. (18) 
p=l 

We denote the probability generating function for the isotropic emission as 

I(t) = e x isotrt-x iBO u _ ( 19 ) 

The total generating function for the probability distribution of photons in pixel is the 
product of the three components 

oo 

J2Pkt k = P(t)-G(t)-I(t). (20) 

k=0 
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The parameters of the point-source distribution in Equation (JSJ) and the level of the isotropic flux 
Xi S otr are found from fitting the probability distribution given in Equation (1201) to the observed 
probability distribution. The details of the fitting algorithm are presented in the next subsection. 

This simple factorized form of the generating function is valid in every pixel. In general, this 
factorization is not possible if we take an average over pixels, unless at most one sources is non- 
isotropic and the other sources can be taken out of the average over pixels. The factorization 
leads to a significant simplification of calculations at the expense that, due to a variation in 
exposure, a part of isotropic EGB may be misinterpreted as a non-isotropic Galactic foreground. 
his effect is expecte d to be small, provided that the Fermi -LAT exposure is sufficiently uniform 



(lAtwood et al.l 120091 ). Another possible source of systematic uncertainty is that the distribution 



of galaxies in the local neighborhood of the Milky Way is non-uniform and a part of emission 
from these galaxies can be misinterpreted as Galactic emission. 



2.4. Fitting Algorithm 

The algorithm has two main parts: 

1. Determine the non-isotropic part Xq &1 of Galactic diffuse emission. We mask the bright 
point sources and use low multipoles of spherical harmonics to find the variable part of the 
remaining flux (Appendix |B]) . Then we add a constant such that min(xQ al ) = 0. In the 
following x^ al is fixed, i.e. we do not vary this component in fitting the pixel counts. 

2. Use pixel counts to find Xi SOtr and dN/dS. The source-count distribution has four parame- 
ters: normalization, break, and two indices. Together with the level of isotropic contribu- 
tion, x iso tr, this gives five fitting parameters which we denote by a = (A, Sb re ak, ni,n 2 , x isotr ), 
where A denotes the normalization in the source-count distribution in Equation ([5]). 

The model prediction for the pixel counts z/ fc (a) is determined by multiplying the right-hand side 
of Equation (J20D by A pix 

A pix • P(t) ■ G(t) ■ I(t) = »k(a)t k . (21) 

k=0 

The expected pixel counts v^{a) are found by expanding the left-hand side of this equation in 
powers of t and picking the coefficient in front of t k . Given the observed number of pixels n k 
with k photons, the likelihood of z^(a) is estimated by the Poisson probability 

n k 

L k = JL-e-"". (22) 



-10- 



The overall likelihood of the model is estimated as 

L(a) = TT^l!! e -^W. (23 ) 

Here the product is over the number of photons k, while the data are represented by the number 
of pixels rifc that contain k photons. This is different from a usual coordinate space fit of maps, 
where the product is over pixels and the data are the number of photons k p in a pixel p. 

We use this likelihood function to find the best-fit parameters a* for a given set of observed 
pixel counts n&. The significance of deviation of a from a* can be estimated as 

^2 = lnL(a,)-lnL(a). (24) 

In the case of large counts the likelihood is well approximated by the Gaussian distribution and 
a (a) is the deviation in sigma values. 



3. Data analysis 

We consider 11 months of Fermi data (August 4, 2008 - July 4, 2009) that were used for the 
first Fermi p oint-source catalog (lAbdo et al.ll2010al ) and for the Fermi analysis of the source-count 
distribution ( lAbdo et al.ll2010d ). 



In order to reduce the PSF, we use front-conve rted gamma-rays, i. e., the gamma-rays con- 
verted in the front part of the Feimi-LAT tracker (lAtwood et al.l 120091 ). with energies between 
1 GeV and 300 GeV. The corresponding (qua dratic) average of the PSF is 0.4° §| For pixelation 
of data, we use HEALPix (IGorski et al.l 120051 ) with the pixelation parameter nside = 32, which 
corresponds to pixel size about 2 degrees J§ We mask the galactic plane within 30 degrees in 
latitude. The number of pixels outside of the mask is N p - lx = 6178 and the number of photons is 
152,143. 

The distribution of pixel counts is presented in Figure [2] on the left. We model this distri- 
bution by a combination of three components: non-isotropic Galactic diffuse emission derived in 
Appendix [Bj isotropic distribution of photons, and the distribution of AGN-like point sources 
described in Section |2~U1 The best-fit models for the isotropic flux and the AGN-like point sources 
are presented in Figure [2] on the right. 



http : //www-glast . slac . Stanford . edu/sof tware/IS/glast_lat_perf ormance .htm 



3 We consider the cases of nside = 16 and nside = 64 in Appendix [Cl 
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Fig. 1. — Average distribution of flux among pixels from Gaussian sources at random positions. The 
width of the Gaussian corresponds to the average point spread function for front-con verted gamma-rays 
with energies between 1 GeV and 300 GeV. The pixels are determined in HEALPix (jGorski et al.ll200a ) 
with nside = 32. x-axis: fraction of the flux, y-axis: density of pixels as a function of / (cf., Equation 
([5])) multiplied by the flux fraction. The fractions are computed by taking 50,000 Gaussian distributions 
at random positions on the sphere. The area under the curve corresponds to the total flux equal to 1. 



In fitting we consider pixel counts with k < 500, i.e., the number of data points is 
-^points = 500. We use £ max = 20 for the Galactic diffuse emission model. The sampling is 
performed with Metropolis-Hastings Markov Chain Monte Carlo method with 1500 steps. The 
average model is presented in Figure |2] and compared with Fermi models in Table [TJ We find 
that the position of the break is almost a flat direction. The results of fitting with a set of fixed 
break positions are presented in Figure [3J 

We find that the contribution of point sources to the gamma-ray flux is (17 ± 2)%. This 
fraction is obtained by integrating the best-fit source-count model presented in Figure [2] and 
in Tabled] from zero to infinity This fraction is not model independent, i.e., it may change for 
different shapes of the source-count function used in fitting, but, provided that the broken power- 
law source-counts distribution gives a good fit to the data, we do not expect a large systematic 
uncertainty due to the change in the shape of the fitting function. In part this can be justified 
by looking at the models with fixed positions of the break: most of the models in Figure [3] have 
the point-source fraction q ps < 0.25. 
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Fig. 2. — Left plot: n/% is the number of pixels with k photons, the red dots correspond to pixel 
counts derived from Fermi data, the errors bars are equal to y / n&. We consider three sources: AGN-like 
point sources (blue dotted line), isotropic Poisson contribution (brown dashed line), and non-isotropic 
Galactic diffuse emission (black dash-dotted line). Note, that the total model on this plot is not the sum 
of the components, the corresponding generating function of the PDF is a product of the generating 
functions in Equation ([20]) . -/Vp i nt s corresponds to the number of points on the x-axis that we have used 
for fitting, k < 500. Right plot: green solid line is the physical model for the AGN-like source counts, 
blue dotted line is the source counts in pixels in the presence of PSF (Equation (fT2l) ). points represent 
expected number of m-photon sources per deg 2 as defined in Equation (jlip . blue star is the value of the 
isotropic flux. 



At smaller position of the break the contribution of point sources becomes larger (g ps ~ 0.3 
for the left point in Figure [3]) but the fit has smaller likelihood. It would be interesting to estimate 
constraints on models invo 
(e.g., 



Abazajian et al. 



2010; 



ving high redshift AGNs t o explain the unresolved part of the EGB 



Neronov fe Semikozl 120111 ) using the statistics of photon counts. 



We can also put an approximate lower bound of 51% on Galactic diffuse emission, and 
an approximate upper bound of 32% on isotropic emission consistent with Poisson statistics. 
In terms of the absolute flux values, the total gamma-ray flux above 1 GeV for |6| > 30° is 
F tot = 1-75 x 10~ 6 s _1 cm _2 sr _1 , the diffuse Galactic flux is F Ga i > 8.8 x 10~ 7 s _1 cm _2 sr _1 , 
the isotropic component of the flux is -Fi SO tr — (5.6 ± 0.6) x 10 -7 s _1 cm _2 sr _1 , the flux from 
AGN-like point sources with dN/dS parameters given in Table Q] and in Figure [2] is Fps = 
(3.0 ±0.4) x 10- 7 s- 1 cm- 2 sr- 1 . 



Using this value of the total flux from AGN-like point sources, we can estimate the gamma- 



10" 9 
■Sbreak (10 



10"' 



Fig. 3. — Results of fitting the AGN-like source-count distribution for fixed positions of the break. The 
significance of deviation is determined from Equation (|24j) as a = \/2AlnL. n\ and 112 are the indices 
above and below the break for differential source-count distribution, q ps is the fractional contribution 
from AGN-like point sources, qi so tv is the fractional contribution from an isotropic source (it provides 
an upper limit on the extra-galactic isotropic contribution additional to AGN-like point sources). Red 
dots correspond to the model used in Figure [2j Green dashed lines provide reference values: n\ = 2.5, 
n 2 = 1.5, q ps = 0.2, and q isotl = 0.3. 



ray flux from undetected sources. Th e total flux in the energy range 1 - 100 GeV from point 
sources at |6| > 30° detected by Fermi (lAbdo et al.ll2010al ) is Fpermi agn — 2.1 x 10~ 7 s"~ 1 cni -2 sr -1 . 
Consequently, the flux from undetected AGN-like point sources above 30° can be estimated as 
^unresAGN = 0.9 x 10~ 7 s -1 cm~ 2 sr -1 . This value is reasonable, since from FigureHJit follows that 
the detected point sources give a good approximation to the best-fit distribution of sources down 
to fluxes below the break in N(S), while the contribution of point sources to the gamma-ray 
background is saturated aro und the break. The extragalactic diffuse emission can be estimated 

do et al.l (hoiObh . In this model, the EGB flux above 1 GeV is 
4 x 10~ 7 s _1 cm~ 2 sr _1 . The contribution of unresolved AGN's to this EGB flux is about 



from the model presented in 
-Pegb 

23%. This value is consistent with the estimations in 



that our analysis is different from 



Abdo et al 



AbdoetaL 



(j2010d j . We would like to stress 



(l2010d j. We do not count any sources, but rather 
do a "blind" statistical analysis. The consistency of results may serve as a non-trivial check of 
both methods. 



-14- 



Table 1: Comparison with Fermi Results (|Abdo et al.ll2010d ). 



Analysis 


ni 


n 2 




CO 

break 


This work 


2.31 ±0.25 


1.54 ± 


0.16 


2.9 ±0.9 


Fermi > 100 MeV 


2.44 ±0.11 


1.58 ± 


0.08 


2.8±0.5 b 


Fermi 1 - 10 GcV 


oq+0.15 
^••3»_0.14 


1.521 


0.8 
1.1 


2.3 ±0.6 



a In units of 10 9 ph s 'cm 2 . 

b Rescaled from the break in lAbdo et al. ( 2010c ) according to energy spectrum ~ E 



4. Discussion 

One of the most interesting problems in gamma-ray astrophysics is to understand the origin 
of EGB flux, which is defined as the total gam ma-ray flux minus the contribution from resolved 



point sources, minus the Galactic foreground (lAbdo et al.ll2010bl ). In general, we can separate 



the sources of diffuse EGB into two big classes: non-Poisson sources and Poisson-like sources. In 



photons during the obser v ation time, these are the AGN-like point sources (e.g., 



1993 



Chiang et al. 



1995 



Stecker fc Salamon 



1996 



Abazajian et al. 



2010 



(e.g.. 


Padovani et al. 


Neronov & Semikoz 



2011I ). In this case the statistics of photon counts in pixels across the sky will be different from 
Poisson statistics due to correlation among photons coming from the same point source. In the 
second case, there is a large number of sources emitting on average muc h fewer than one photon 



during the observation time, e.g., star- fo rming and star-burst galaxies (jPavlidou fc Fields 



Thompson et al 



2007 



Fields et al 



2002 



2010j). The statistics of photons in this case will be very close 



to Poisson statistics. 

These two possi bilities are relatively hard to separate based only on energy spectrum argu 
ments (see, however, 



Abazajian et al. 



2010 



Stecker fc Venters 



201 



0;' 



Neronov fc Semikoz 



201l|). 



Alternatively, on e can use a correlation between gamma-ray emission and multi- wavelength emis 



Li&Cao 



Hop kins et al 



2007 



sion from AGNs (lUrrv fc Padovani 
2010l: 



1995 



Miicke k Pohl 



2000 



Abdo et al. 



2010d; 



Angela kis et al. 



2011) together with AGN population studie s fe.g.. iPadoyani fc Giommi 



Padovani et al. 



2007 



Rigby et al 



2008 



Hodge et al. 



2009 



1995 



Smolcic et al. 



20091 ) to infer the AGN contribution to gamma-ray background. 



In this paper we show that the statistics of photons in pixels can also be used to constrain the 
distribution of point sources below the detection limit. We find that there is a slight preference 
for two populations of point sources: the AGN-like point sources with a break at relatively 
high flux and a population of faint sources such as star-forming galaxies rather than a single 
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Large scale model variation 
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Fig. 4. — Left plot: different lines correspond to different models for Galactic diffuse emission of 
gamma-rays corresponding to £ max = 10, 20, 30 (Appendix [B]) . Solid green line corresponds to the 
range of photon counts (1 to 500) used in fitting. Dashed green line is the extrapolation of the model. 
Magenta triangles c orrespond to total source counts found by Fermi collaboration above 100 MeV (right 
plot in Figure 9 of 



AbdoetaL 



2010d ) with flux rescaled assuming energy density spectrum ~ E~ 2A . 
Blue circles correspo nd to sources found by the Fermi collaboration in 1 — 10 GeV energy bin (center 
plot in Figure 17 of lAbdo et al.l ( 2010c)). Green squares correspond to 1 — 100 GeV flux for sources 



in the first Fermi catalog (jAbdo et al.ll2010al ). Right plot: models with different fixed positions of the 
break presented in Figure 



population of AGN-like point sources with a break at a smaller flux. This statement is insensitive 
to the energy spectrum since we use the integrated flux above 1 GeV but it may depend on the 
form of the function used to fit the point source number counts. Possible sources of systematic 
uncertainty include the source smearing due to PSF, non-isotropic Galactic diffuse emission, non- 
isotropic emission (on large scales) from nearby galaxies, and non-homogeneous exposure. The 
interpretation of the point-source number counts may also be affected by clustering of sources 
on scales smaller than the PSF. 

A stronger constrain on the population of gamma-ray point sources may come from an 
unpixelized analysis of the data. In this more general analysis, the likelihood of a model de- 
pends on the full gamma-ray data (positions on the sky, energy, arrival time) rather than the 
counts of gamma-rays in pixels. This likelihood would have the form of an integral over flux 
times the source number counts, times an integral over all possible positions of the sources on 
the sky. The photon counts approximation is simpler computationally and already gives rough 
constraints on the source population, while the full unbinned analysis may provide stronger con- 
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straints at the expense of computational complexity. Examples of unbinned analysis using some 
part of gamma-ray data include the gamma -ray two-point correlation function 



2009 



Geringer-Sameth &: Koushiappasll2010l ). a ngular power spectrum (lAndo fc Komatsu 



Hensley et al.ll2010l ). nearest neighbor statistics (jSoltanll201ll ). etc. 



function (e.g., 


Ave et al. 


(Ando & Komatsu 


2006; 



We also believe that the techniques of generating functions in the study of photon counts 
statistics, known for some time in radio and X-ray observations, may have further important 
applications in data analysis of current and future radio, IR, X-ray, and gamma-ray observations. 
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A. Derivation of statistics 



In this Appendix we derive the photon counts statistics in the presence of point sources. The 
derivation is based on the following well known property of generating functions for probabilities: 
the generating function of a sum of two ind ependent sources is the product of the corresponding 
generating functions (e.g., iHoel et al.lll97ll Section 3.6). Indeed, consider a box where we can 
put red balls with probabilities p^, k — 1, 2, 3 . . ., and blue balls with corresponding probabilities 
qk, k — 1, 2, 3..., then the probability to find k balls of any color is 



Pk>qk~k', 

fe'=0 



k = 1, 2, 3... 



(Al) 



The last equation is the rule for multiplication of polynomials or power series. Let 



(A2) 



then 



R(t) = P(t)-Q(t). 



(A3) 



-17- 



As before, denote the probability to observe k photons in a pixel by p^. The corresponding 
generating function is 

P{t) = Y,Pkt k . (A4) 

k 

The knowledge of the generating function is equivalent to the knowledge of all pk, provided that 
every pk can be found from P(t) by picking the term in front of t k . 

We will now derive the generating function for probabilities to observe k photons from a 
collection of point sources. Denote by x m the average number of point sources inside a pixel that 
emit exactly m photons during the time of observation. We assume that the probability to find 
n m m-photon sources in a pixel is given by the Poisson statistics 

PnJ*™) = -V"*" , (A5) 
I'm- 

then the probability to find k photons from m-photon sources is 

p ( m ) _ I Pn m ( x m), if k = m ■ n m for some n m ; 
^ k 1 0, otherwise. 

The corresponding probability generating function is 

p (m) ^ = Pk t k = tm ' nm ^je~ Xm = e Xmtm ~ Xm . (A7) 



The generating function of probabilities pk to observe k photons from any sources is the product 
of the generating functions for every m 



oo 



E^* 1 = n 

^ ] (x m t m — x m ) J . (AJ 



e 

A-i") m=l 

oo 

exp 



\m=l 



If we denote t = e and 



P(cj) = X>e 2 ™ fc , (A9) 

fc=0 

oo 

X{u) = 5> m e 2 ™, (A10) 

m=0 

then Equation flA8j) can be represented as 

P( w ) = e*< w >-*<°>, (All) 
which is the d iscrete Fourier transform analog of characteristic functions considered in, e.g., 



Scheuerl (Il957l ). 
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Angular power spectrum, E = 1 - 300 GeV 

10 4 F 1 1 1 1— " 



C5* 




10° 10 1 10 2 10 3 

I 

Fig. 5. — Blue line (upper): angular power spectrum of gamma-ray data with |6| > 30°. Red line 
(lower): same as the blue line, but with additional masking of Fermi gamma-ray point sources (Figure 
E]). The normalization is chosen such that for the Poisson noise (Cj) = 1 (constant green line). C/'s 
below i ~ 10 are dominated by the variation in Galactic diffuse emission on large scales. Constant C/'s 
above t = 10 are due to contribution of point sources. Decay of Cj's above I ~ 100 for the blue line is 
due to detector PSF. In our case (PSF) « 0.4°, which corresponds to I > 400. A smoothed model for 
large-scale structure distribution is obtained by spherical harmonics decomposition up to I = 20, the 
corresponding map of the model is presented in Figure El 



B. Non-isotropic Galactic diffuse emission 

In this appendix, we describe a model for the non-isotropic Galactic diffuse emission. This 
signal is "filtered" by low multipoles of the gamma-ray data. In this approach, the homoge- 
neous part of Galactic emission is indistinguishable from extragalactic flux. Also some part of 
extragalactic emission from galaxies and galaxy clusters close to Milky Way may lead to a signal 
that varies over the sky and can be misinterpreted as non-isotropic Galactic emission. How- 
ever, the majority of extragalactic emission comes from higher redshifts where the distribution 
of extragalactic sources is sufficiently uniform on large scales. 

In Figure El we plot the angular power spectrum for the data at high latitudes (161 > 30°) 
before and after masking the point sources in the first Fermi catalog (lAbdo et al.l l2010af ) . The 
angular power spectrum is 

^airrEi^i 2 ' (B1) 
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Fig. 6. — Sky maps for photon counts and for a Galactic diffuse emission model derived from the 
photon counts by using the spherical harmonics decomposition with I < 20 (Equation (IB2I) ). Residuals 

i I model— data! 

are m sigma values: J — , ■ , — '-. 

V model 



where dim s circ the spherical harmonics coefficients of f(p) ■ m(p), where f(p) is the number of 
photon counts inside pixel p and m(p) is the mask function. The mask function is equal to one 
(zero) for \b\ > 30° < 30°). In the case of masked point sources, we also have m(p) = when 
a Fermi point source is inside the pixel or within two PSF from the boundary of the pixel. We 
also subtract the average of the data within the unmasked region in order to avoid a non-trivial 
contribution from a constant source inside the window. We choose the normalization such that 
for the Poisson noise (Cj) = 1. 

The algorithm for estimating the component that varies on large scales is as follows: 
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k-photon pixel counts, no PS 



Pixel counts, no PS 
Diffuse model 




Fig. 7. — Photon counts vs smooth model with Fermi point sources subtracted. The Galactic diffuse 
emission model is obtained by taking the spherical harmonics of the data with i < 20. The model 
photon counts are derived from the probability generating function in Equation (|16|) with x^ie given in 
Equation (lB2l) . 



1. Calculate the aim's inside a window larger than needed (in order to avoid edge effects). We 
use |6| > 20° for the data above 30°. We also fill in the pixels with point sources by the 
average of nearest neighbors. 

2. The large-scale distribution of gamma-rays is defined as 

^liiax ' 

Xdiff(p) = A Yl a l™Ylm{p) + B, (B2) 

1=0 m=-l 

where we find the coefficients A and B from the best fit of the diffuse model to the photon 
counts in the pixels without point sources (see Figure [7]) . 



3. We represent B = B min + x isotr so that 



( p ) = AJ2J2 a imYim(p) + B 



(B3) 



1=0 m= 



is non- negative for |6| > 30°. In fitting to full data, A and B inin are fixed, while a; isotr is 
allowed to vary together with parameters describing the AGN-like point sources. 

An example of the diffuse emission model £diff(p) for Z max = 20 is presented in the middle plot of 
Figure O The top plot represents the counts of photons in pixels inside the window |6| > 30° with 
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5 >1.0GeV (P hs 2 ) 



Fig. 8. — Same as in Figure [2] but with the pixel size of about 1° corresponding to HEALPix parameter 
nside = 64. 



masked Fermi gamma-ray point sources. The bottom plot represents the deviation of the model 
from the data in "sigma" values. In Figure 0] we study the effect of changing / max = 10,20,30. 
The difference is rather small, i.e., already Z min = 10 captures the large-scale distribution of 
gamma-rays reasonably well. 



C. Variation of pixel size 

In the data analysis in Section |3l we use the pixel size of about 2° corresponding to HEALPix 
parameter nside = 32. In this appendix, we repeat the analysis of Section [3] for different sizes of 
pixels, in particular we consider nside = 64 and nside = 16 corresponding to pixel sizes 1° and 
4° respectively. 

The results for nside = 64 are presented in Figure [HJ The total number of pixels inside 
the window below and above 30° in latitude is about 25,000, i.e., the statistics of pixel counts 
is relatively good, but the effect of PSF becomes more significant than in the case of nside = 
32. For PSF = 0.4°, less than about 60% of the flux from a source can be inside one pixel. In 
particular, one can note that the number of sources with small fluxes is about two times larger 
than in the input model. Nonetheless, in the case of nside = 64 the best-fit model is similar to 
nside = 32 case. The position of the break, Screak — 1-1 x 10~ 9 ph s _1 cm~ 2 , is somewhat lower 
than the values in Table Q] but the overall source count distribution (Figure [9]) is consistent with 



the source counts found by Fermi collaboration (lAbdo et al. 



2010a. 
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Fig. 9. — Comparison with Fermi source counts (jAbdo et al.l l2010al Jd) in the case of analysis with 
HEALPix parameter nside = 64 (Figure [8|) . 



The results of fitting for nside = 16 are presented in Figure [TOj In this case, the PSF can 
be neglected. However, the number of pixels inside the window is only about 1,500. As one can 
see on the left side of Figure [TTJ], the error bars are rather large and the photon counts in most 
of the pixels are dominated by the contribution from diffuse emission. As a result, our method 
is not sensitive to the contribution from point sources for nside = 16. In particular, the MCMC 
is dominated by models with large Screak- The best fit value S^eak = 1-3 x 10~ 8 ph s _1 cm~ 2 is an 
order of magnitude larger than the values in Table [U The reason is that models with large 5b re ak 
have a comparable likelihood to the models with small Sbreak but there are many more ways to 
choose a large Screak- 

We conclude that our method is relatively robust to the change of pixel size: the cases of 
nside = 32 and nside = 64 give very similar results, unless the number of pixels is not sufficient 
to draw any conclusion about the statistics of sources, as in the case of nside = 16. 
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n x =5.83 ±3.15 




n 2 =1.60 ±0.07 







Fig. 10- 
nside = 16 
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Same as in Figure[2]but with the pixel size of about 4° corresponding to HEALPix parameter 
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